"nevtot: total number of events for all individuals"
CALC nevtot = nval(indiv)
CALC nevtot2 = nevtot - 1
"nindivs: number of individuals"
CALC nindivs = nlevels(indiv)
"lex: number of risk group cut points"
CALC lex = nval(ex)
"lhist: length of history variate (all cut points for 1 individual)"
CALC lhist = lex + nval(t) + 2
"lvars: number of intervals for 1 individual (inc. 0 intervals)"
CALC lvars = lhist - 1
"lall: number of intervals for all individuals"
CALC lall = nindivs*lvars

"z: non-zero if the first event for each individual"
VARI [nvalues = nevtot] z 
CALC z$[1] = 1
CALC z$[2...nevtot] = indiv$[2...nevtot] - indiv$[1...nevtot2]

"restrict sta, end, ex to the first record only"
SUBSET [condition = z .NE. 0] oldvector=sta,end,ex[]

"if covariates exist, restrict to first record and set up variates 'cov'"   
CALC ncovs = nval(covariates)
IF ncovs .GT. 0
SUBSET [condition = z .NE. 0] oldvector=covariates[]
VARI [nvalues = lall] cov[1...ncovs]
ENDIF

"set up variates for the following:\ 
individual: factor for each individual\
exposure_group: exposure group risk factor\ 
age_group: age group factor\ 
interval: interval length\ 
nevents: no. of events\
(1st 3 will be converted to factors before model fitting)"
VARI [nvalues = lall] individual,exposure_group,age_group,interval,nevents

DELETE [redefine = yes] nevtot,nevtot2,z,lall

"loop for each individual" 
FOR ind = 1...nindivs; 

	"pos1 & pos2: where to place data for each individual"
	CALC pos1 = ind * lvars - lvars + 1
	CALC pos2 = ind * lvars 

	"factor for each individual"
	CALC individual$[pos1...pos2] = ind

	"covariates - if they exist"
	IF ncovs .GT. 0
	FOR i = 1...ncovs
	CALC cov[i]$[pos1...pos2] = covariates[i]$[ind]
	ENDFOR
	ENDIF

	"exposure risk group cut points"
	VARI [nvalues = lex] excuts
	CALC excuts$[1...lex] = ex[]$[ind]

	"hist: ordered cut points (start, end, age and risk group)"
	POINTER [VALUES=(sta, end)$[ind], excuts, t] histp
	VARI [nvalues = lhist] hist
	EQUATE histp; newstructures = hist
	SORT hist

	"number of events in each interval" 
	SUBSET [condition = indiv .EQ. ind] oldvector = eventday;\ 
		newvector = eventdayi
	TALLY [print=*] eventdayi; limits = hist+0.5; frequencies = neventsi
	CALC nevents$[pos1...pos2] = neventsi$[2...lhist]

	"exposure group risk factor"
	GROUPS [lmethod = give] vector = hist; factor = exi;\ 
		limits = excuts+0.5; levels = exc
	CALC exposure_group$[pos1...pos2] = exi$[2...lhist]

	"age group factor"	
	GROUPS [lmethod = *] vector = hist; factor = agegri; limits = t+0.5
	CALC age_group$[pos1...pos2] = agegri$[2...lhist]

	"interval lengths"  
	CALC interval$[pos1...pos2] = hist$[2...lhist] - hist$[1...lvars]
	"set interval = 0 if cut points are less than sta or greater\
	than end (these will be taken out before fitting the model)."
	FOR i=1...lvars; k=pos1...pos2
		IF hist$[i] .LT. sta$[ind]
		CALC interval$[k] = 0
		ENDIF
	ENDFOR
	FOR j=2...lhist; k=pos1...pos2
		IF hist$[j] .GT. end$[ind]
		CALC interval$[k] = 0
		ENDIF
	ENDFOR

ENDFOR

DELETE [redefine = yes] covariates,ex,lex,histp,eventday,eventdayi,indiv,\
neventsi,exi,excuts,exc,agegri,t,hist,sta,end,pos1,pos2,lhist,lvars

"cut out all intervals in which the interval length is 0 and\ 
convert covariates to factors"
IF ncovs .GT. 0
RESTRICT individual,interval,age_group,exposure_group,nevents,cov[];\ 
		condition = interval .NE. 0
GROUPS [redefine=yes] cov[]
ELSE
RESTRICT individual,interval,age_group,exposure_group,nevents;\
		condition = interval .NE. 0
ENDIF

"log the intervals"
CALC loginterval = LOG(interval)

"convert individual, agegr and vacc to factors"
GROUPS [redefine=yes] individual,age_group,exposure_group

"... and finally fit the model :-)"
MODEL [distribution = poisson; link = logarithm;\
	 offset = loginterval; groups = individual] nevents
FIT [print = model, summary] ##fit
RKEEP Y = nevents; estimates = estimate; se = s_e; lower = CI_low;\
	upper = CI_hi
CALC RR = exp(estimate)
CALC CI_low = exp(CI_low)
CALC CI_hi = exp(CI_hi)

"print out estimates, SEs, RRs and CIs for RRs"
PRINT estimate, s_e, RR, CI_low, CI_hi; decimals = 3

"table of total number of events by age group and vaccine risk factor"
TABULATE [print = totals; classification = age_group,exposure_group] nevents